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Introduction 

In the present paper I will describe how the first variant of the "Godunov's scheme" has 
been elaborated in 1953-1954 and tell about all modifications realized by myself (until 
1969) and the group of scientists from the Institute of Applied Mathematics in Moscow 
(which has become the M.V.Keldysh Institute of AppUed Mathematics). 

At the time these modifications (see Sections 2,3) were carried out, other algorithms 
were developed, in particular second order schemes for gas dynamics problems with a 
small number of strong and weak discontinuities [1-3]. We performed many calculations 
based on the first codes written by V. V. Lucikovich. More complicated problems resulted 
in an elaboration of very artful approaches to divide the whole computational domain into 
sub-domains which have been developed by A.V.Zabrodin. This procedure resulted in the 
necessity to develop algorithms of grid construction. 

In 1961-1968 G. P. Prokopov and I carried out approaches to the construction of moving 
grids which were used in serial calculations by A. V. Zabrodin, G. N. Novozhilova and 
G. B. Alalikin (see [4-6]). The problems appearing in the grid construction forced us to 
solve elliptical systems (see [7-8]). The methods elaborated here, have later been employed 
in elliptical spectral problems and have been presented in my papers on numerical linear 
algebra (see [9,10]). The number of interesting observations made during the analysis 
of my calculations gave many discussions at the Moscow University and, after 1969 - 
at the Novosibirsk University. As a result of such discussions the criterion of spectral 
dichotomy [11,12] has been developed and high-precision algorithms to calculate singular 
vectors have been constructed (see [11,13]). It is difficult to imagine that the reason of 
such investigations has been the elaboration of approaches to the gas dynamics calculus 
and numerical grid constructions. 

During my studies at Moscow University I learned the differential equations theory in 
seminars of I. M. Gelfand and I. G. Petrovskii. The latter focused my attention on gas 
dynamics problems and proposed to me to use stationarization methods to study transi- 
tional flows (with sub- and supersonic regions). My qualification work was devoted to the 
stationarization of a flow inside a nozzle (however, only in subsonic regime and with artifi- 
cially added time derivatives introduced in Chaplygin's equation). Petrovskii's idea about 
stationarization in practical form was published in 1961 (see [6]). The technical statement 
was presented in the qualification work of G. P. Prokopov performed under my supervision. 
The coding was made by G. N. Novozhilova. 

The elaboration of numerical schemes was carried out at the same time with attempts to 
have a better understanding of the notion of generalized solutions to quasi-linear systems 
of equations. As a rule, the hypothesis about possible definitions and properties of the 
generalized solutions were preceded to the construction of numerical schemes, which used 
these properties. At the same time, I tried to prove the formulated hypothesis. To my 
deep disappointment, these attempts had no success. On the contrary, they often led to 
contradictory examples. But, at the same time, a numerical scheme, more precisely its 
modification, based on using the Euler coordinates, moving grids and tracking methods for 
strong and weak discontinuous, was used in customary calculations. 

1 How the scheme has been elaborated 

In autumn 1953 M. V. Keldysh and I. M. Gelfand proposed me to elaborate a variant of a 
method suggested by J. Von Neumann and R. D. Richtmyer It consisted in the introduction 
of artificial viscosity in the gas dynamics equations. The goal was to finish the algorithm 
for serial engineering calculations by the spring of 1954. By this time, the first electronic 
numerical device "Strela" had to be delivered at our Institute. 

I knew that one variant of the required algorithm in our Institute had been prepared by 
A. I. Zhukov. I had the opportunity to study it and the goal was to modify this algorithm 
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or to suggest a new variant. From the point of view of the administration of our Institute 
it was necessary to continue the study in this direction to be able to solve efficiently this 
problem to a fixed date. 

The method of A. I. Zhukov exactly coincided with the scheme published by P. D. 
Lax one year later [14]. The paper I received from A. I. Zhukov contained some ideas 
and hypotheses on which the numerical scheme was based. Some of the hypotheses were 
erroneous but during my first reading of the paper I was confident in them and decided to 
rely on these hypotheses in my work. 

The principal difference between Zhukov's method and the approach by J. Von Neu- 
mann and Richtmyer lies in the use or not of the artificial viscosity. A. I. Zhukov under- 
stood during his experiments that smearing of shock waves was made automatically be- 
cause of discrepancy between approximate and exact differential gas dynamics equations. 
A. I. Zhukov explained this discrepancy by analyzing the truncation terms of the difference 
scheme, which present numerical viscosity, and formulated the idea of the first differential 
approximation. 

He began to use the second order difference schemes but faced with strong oscilla- 
tions near the shock wave front. To explain this fact he used exact solutions based on Airy 
functions of linearized equations in which numerical viscosity was introduced to model 
truncation terms of the numerical scheme. The experiments resulted in the conclusion 
that the second order schemes had to be avoided, and A. I. Zhukov concentrated on the 
development of the first schemes for the simplest quasi-linear Burgers equation and the 
gas dynamics equations in the Lagrangian coordinates. The whole attention was concen- 
trated on the investigation of numerical solutions of problems with a steadily moving shock 
front. When the first differential approximation was used, the problem was reduced to the 
solution of the ordinary differential equations which could be compared with numerical so- 
lutions. Based on these comparisons A. I. Zhukov formulated a necessary condition which 
any difference scheme should satisfy. This condition means that the first differential ap- 
proximation should have the form of conservation laws, however it was formulated as a 
hypothesis of a necessary and sufficient condition. 

I did not doubt in that hypothesis and put it in the cornerstone of my attempts to build 
a scheme without oscillations but with second order accuracy. To reach this goal I began to 
compare different second order schemes which are not conservative but they are conserva- 
tive at the first order. The comparison of different variants was based on the same solutions 
which were important for a whole set of typical problems. As a result a rather satisfactory 
variant was chosen. Three months of intensive work were spent for such a choice. N. M. 
Zueva helped to perform calculus, and numerical experiments were performed by V. V. 
Paleichik and two of her students who studied technology of calculation based on mechan- 
ical adding machine "Mersedes". M. V. Keldysh and I. M. Gelfand were interested in these 
investigations; they listened to reports about the current state of this work biweekly. 

The next stage of the investigations was to prove the efficiency of the scheme based 
on various numerical experiments with different spatial steps and different solutions. But 
already during investigation of the dependence of the solution on the spatial step we faced 
with the situation that completely crossed out all results of our three-month work. 

To clarify the difficulties we have faced I have to teU about one hypothesis from Zhukov's 
report. He compared the solution of a problem about steadily moving shock waves by us- 
ing the first differential approximation with the real profile obtained from the numerical 
calculations. They were different and he supposed that this difference would be decreased 
if the spatial step tends to zero. 

I decided to check this hypothesis and asked V. V. Paleichik to make corresponding 
calculations. One set of calculations with one step size was made by the first student, and 
the calculations with one-half of that step by another. Rapidly, in five-ten minutes V. V. 
Paleichik reported me about senseless of the proposed experiment. The results of both 
experiments were absolutely the same, and obviously a different result was not possible 
because spatial step itself did not play any role in calculations. Only the relation between 
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time and spatial steps took part in the experiments and it was chosen constant to keep 
Courant's number constant. So, finite-difference profiles were not modified in the process 
of step decreasing and did not converge to the dimensionless profile obtained for the first 
differential approximation. At the same time, the comparison of the velocity of the exact 

and finite-difference shock fronts demonstrated no significant difference (« 3 4%), but 

this difference did not decrease during steps reduction. 

Obviously, this situation put me under stress and forced me to change my previous point 
of view. 

Firstly, it was obvious that the numerical scheme should guarantee the correct velocity 
of a steadily moving shock wave and correct values on both sides of the front even if the 
step size equals to 1. The hope, that the accuracy in calculation would be high if the step 
was small enough, disappeared. 

The conclusion was to use exact conservation laws without simplifications related to 
the use of conservation laws obtained from the first differential approximation. For this, 
the approximation of the equations 

du dp iv ) 

1 — = , 

dt dx ' 

dv du p 
dt dx ' 

should be done to provide the following conservative form of difference equations 

- =0, 



At Ax 

Jjrn _ jjn 
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To describe the velocity field one should use both and U^^i- Similarly, the specific 
volume V and pressure p(v) related to it are approximated by both v'^ and P™ i respec- 
tively. The different formulas connecting the lower and upper case letters result in the 
different schemes. The resulting scheme should provide the smoothness of a numerical 
profile. Based on this scheme and numerical conservation laws one can conclude that ve- 
locity of the wave is correct and values on the right and left sides of the shock are connected 
by the known Hugoniot relations. 

To construct a scheme I decided at first to obtain a smooth front in the Umit case of weak 
waves, i.e. in the case of acoustic equations, which can be considered as the gas-dynamics 
equations with the simplest equation of state p{v) = v (a^ = const. > 0). In this case, 
the shock waves are transformed into the discontinuities of the Riemann invariants u±av, 
and the equations are reduced to a canonical form 

d{u±av) d{u±av) ^ 
dt dx 

I decided to choose numerical schemes which allow (at least in this simplest case p(v) = 
—a^v, a = const.) the reduction of these equations to independent equations for the Rie- 
mann invariants u± av. I used the simplest scheme conserving the monotonicity of the 
Riemann invariants which was known and widely used at that time. The monotonicity gua- 
ranteed a required smoothness. At that time 1 did not have another criterion and there was 
no more time for proposing something else. "Strela" was already being installed in our 
Institute. At that time 1 could show for the linear case that only the first order schemes 
conserve the monotonicity. 1 decided not to provide further efforts in order to to construct 
second order schemes. 
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The chosen scheme in the simplest case had the following form : 

T ^ 2h ^ 2h " ' 

T 2/1 2/i 

and transfomed into a finite-difference conservation laws after introduction of the follow- 
ing notations 

1 I — Q — -p d ~ — T ~r d 



"+3 2 2 2 2 

"+i - 2 2 ~ 2 2^^ • 

It was decided to use the scheme in a general nonlinear case, but the constant a in formu- 
las P^_^i/2 - ^'^'^ replaced by some, generally speaking, non-constant values 
^n+i/2- ^^^^ variant to calculate a™_,_i/2 '•^^ following 



\ (V-p'(^^) + V-p'«Vi)) ' 



These formulas were tested in the numerical experiments during which the monoto- 
nicity of the shock profile was checked. The initial data for the shock were chosen as 
constant values on the right and left sides of the discontinuity. These constant values were 
connected by Hugoniot relations. 

During these experiments 1 began using other formulas for 0^^1/2- I" particular, 1 
wanted to conserve monotonicity at the first step and 1 could construct some interpolation 
formula for a™_^i/2 which gave the most acceptable results. 

Here 1 should add some precisions to my story to correct the simplifications 1 intro- 
duced. We soon began the above described experiments not with the simplest gas dynamics 
equations we used before but with the system of three equations 



du dp{v , E) 
dt dx 
dV_du 
dt dx 



= 0, 
= 0, 



with the following form of numerical conservation laws 
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A selection of the interpolation formulas for a™_)_i/2 became in this case a selection of 
variants to calculate ^ function of u™, -y™, E™, 
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The adequate variant has been found after 2-3 weeks. It is important to remark that in the 
chosen formulas the role of aj^_)_j^2 played the expressions 

/ (7 + l)Pn+l/2 + (7-lK 
V '"n 

which depended on the quantity Pn+1/2 (in these experiments the equation of state E = 
{pv)/{"f — 1) was used). It occurred to me that I saw such expressions somewhere in the 
Uterature. 

Not long before to that the book of L. D. Landau and E. M. Lifshitz "The Mechanics 
of Continuous Media" was published. In this book these quantities took part in formu- 
las which described the solution of the Riemann problem. I had to modify slightly my 
interpolation formulas to obtain the final form of my scheme. These modifications were 
significant only in the case of strong rarefaction waves. It happened in March 1954. At 
the end of March, V. V. Lucikovich was included in our group to write a code according to 
technical tasks that we composed by taking into account useful reconnmendations of K. A. 
Semendjaev. At the end of April 1 was on holidays for two weeks, and on the 5th May we 
began our experimental and applied calculations. 

On the 4th of November 1954 I defended my PhD thesis containing the description 
of the suggested scheme. But this work was published only in 1959. Just before, I tried 
to pubUsh the work in some journals without success. The journal "AppUed Mathematics 
and Mechanics" refused to publish it because it was purely mathematical and without any 
relation to mechanics. In one mathematical journal (I do not remember which one) the 
refusal was motivated with the opposite statement. After that I. G. Petrovskii as a member 
of the editorial group helped to publish it in the journal "Mathematicheckii Sbomik" (see 
[15]). 

I read the paper by P. D. Lax [14] only after my PhD thesis defense. The scheme de- 
scribed in this work coincided with Zhukov's scheme but without questionable hypotheses 
which were in Zhukov's report and whose analysis helped me to construct my scheme. If 

1 had received this paper one year earlier, the "Godunov's scheme" would never have been 
constructed. 

2 Problems of approximation and effective accuracy 

Indeed, the work on modifications of the scheme invented in 1954, was actively pursued. 
Consumers were very suspicious toward my scheme. 

In our Institute, more precisely, in Steklov's Institute, from which Keldysh's Institute 
was separated in 1953, the basic hydrodynamic calculations were organized by K. A. Se- 
mendjaev much before I came. They were made with the help of the mechanical adding 
machine "Mersedes" using the method of characteristics by a large group of colleagues. 
The algorithm was thoroughly thought over by K. A. Semendjaev and A. I. Zhukov. The 
basic attention was payed to the table of results and the reduction of intermediate records. 
This reduction was provided with memory cells of the "Mersedes". All results were pre- 
sented in graphical form and checked very carefully. 

The technique of the calculations was absolutely automatic. My female colleagues 
worked on machines as pianists and simultaneously discussed their household and other 
typical problems of female interest. They mocked especially young researchers among 
whom I was. We were named by the contemptuous word "that science". They supposed 
we built pseudo-scientific theories, useless to overcome difficulties which only led to the 
delay of calculations. One should remark, their salary depended on the calculation volume 
which was defined by the number of table rows without mistakes. 

The graphical presentation of results on a graph paper had a very high quaUty, and 
contained a visual plot of velocity, pressure and density fields, and also the trajectories of 
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contact discontinuities and characteristics giving the possibility to follow the domain of 
influence of the initial data. 

Our customers - physicists and engineers - had a habit to a perfect form of the infor- 
mation obtained and it was a reason of their displeasure toward the first calculation results 
based on numerical schemes. Now it is well-known that the numerical calculations of the 
discontinuous gas dynamics solutions present a "numerical" micro-structure which looks 
like errors providing a caricatural picture of the real flow. The comparison of computer and 
"hand-made" calculations resulted in the opinion that the code had errors or that the initial 
data were erroneous. We had to check carefully all situations, we had to explain the results 
and to modify either the results or the calculation scheme. 

I will describe now two such examples, the most interesting from my point of view. 

During calculations based on my scheme of strong isentropic rarefaction waves one 
could observe a significant growth of the entropy which was the source of distrust toward 
the numerical results. One noticed a substantial reduction of such effect only after 3-4 
years (in 1957-1958) when we studied 2-D problems. Such problems, as a rule, do not 
have a simple formulation in the Lagrangian coordinates and we had to use the Eulerian 
description. Simultaneously, numerical grids became comphcated, they became moving 
ones. The main idea - to use the exact solution of the Riemann problem as an element of a 
numerical scheme - has been applicable in 2-D problems also, but not in a transparent way. 
The first test calculations of 1-D problems earlier computed were repeated using new 2-D 
codes. We were very surprised by the fact that 2-D codes based on the Eulerian coordinates 
resulted in essentially less parasitic increasing of the entropy in the rarefaction waves. 

The reason of the effect is that the Lagrangian coordinates (even in the simplest 1-D 
case) are unfit to the formulation of the generaUzed solutions of the gas dynamics equations. 
These equations admit the formation of vacuum regions which are treated in the Lagrangian 
coordinates by singularities of the delta function type. Usually such singularities are not 
supposed to be present, and they are modelled very badly numerically. The analysis of the 
results showed that even a sufficient decreasing of the step size in Lagrangian coordinate 
in a strong rarefaction zone does not result in a noticeable distance reduction between Eu- 
lerian coordinates of neighboring nodes. I used this fact in 1961 in the theoretical work 
[16] to prove that my scheme approximates (if the Eulerian coordinates are used) the gas 
dynamics equations in a sense of conservation laws (1-D case). At the beginning of this 
work I hoped to prove a convergence to the exact solution but I failed to that. Even the 
order of approximation appeared to be different from one, it was equal to 2/3. Remember 
that, formaUy, my scheme has the first order of approximation based on the first differ- 
ential approximation. I should remark that indeed the order of approximation is probably 
higher than 2/3. Apparently, based on estimates of the variations of solution suggested 
by J. Glimm [17] one can prove that this order is equal to one. Unfortunately I did not 
investigate this issue in details. However at the end of the 50th I was interested basically 
in the accuracy of the approximated solutions, and not in the problem of approximation of 
equations (i.e. what power of the spatial step evaluates the error - the difference between 
the calculation result and the exact solution). The obtained estimates for approximation 
played a preliminary role for theoretical investigation which I could not finish. Despite 
of it, in 1957-1958 V.S.Ryabenkii and I carried out an experimental study of the conver- 
gence of approximate solutions calculated with my scheme to the exact solution. These 
observations about calculation results allowed me to conclude that the convergence had to 
be considered as weak and not strong to exclude the influence of significant deviations in 
one-three nodes. As a rule, the occurring of these deviations is related to the interaction of 
two shock waves or reflection of a shock wave from a contact discontinuity. 

During our discussions, V. S. Ryabenkii suggested an approach which permitted us a 
reduction of the study of weak convergence to the study of strong convergence in 1-D pro- 
blems. For this, one needs to compare not quantities controlled by conservation laws but 
integrals (or even multiply integral) of these quantities with respect of spatial coordinate. 
We supposed to demonstrate the first order of weak convergence, i.e. the order correspond- 
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ing to a formal order of approximation. We were very surprised to see that the observed 
effective convergence order was less than 1 in experiments with centered expansion waves 
or large gradients in smooth domains restricted by strong discontinuities. 

At that time we were very busy and had no motivation to prepare a detailed paper, 
especially because the conclusions were unpleasant to us (to me at least) and did not cor- 
respond to our expectations. 1 only made a brief presentation at a scientific conference at 
the Moscow University. From my point of view, nobody put attention on this presentation 
except N. S. Bachvalov who began the theoretical study of the convergence rate for the 
solutions for the Burgers equation as e ^ 

du du d^u 
dt dx dx^ ' 

and obtained the estimate of the error e In jej. At the same time, all my attempts to attract 
attention of researchers and engineers to these effects had no success. 

My colleagues A. V. Zabrodin, K. A. Bagrinovskii, G. B. Alalikin and I were involved 
at that time in preparing 2-D calculations. That is the reason why the question if the first 
order scheme has the first order convergence rate, or, more precisely, what is this order, was 
not understood. 

I know that later B. van Leer, P. Woodword, P. CoUela, PD. Lax, A. Harten and many 
others (see for instance [18-21]) suggested to construct efficient numerical schemes with 
order of approximation higher than one. At that time I had already left Keldysh Institute in 
Moscow and moved to Novosibirsk. Here I was involved in other problems. Among these 
issues the basic place was devoted to problems which arose during my work in Moscow 
and they required deep theoretical investigations. They were, in particular, issues concern- 
ing an accurate mathematical formulation of the conversation laws and thermodynamical 
identities, energy integrals for the hyperbolic equations, the mathematically correct formu- 
lation of the elasticity theory, and the statement of problems in numerical methods of linear 
algebra. I had already mentioned that this latter issue appeared from attempts to construct 
2-D numerical schemes. The investigation of these various and interesting problems (see 
for instance [22,23]) did not allow me to continue the study of numerical methods and a 
deep understanding of modem numerical schemes for the gas dynamics equations, some 
of which named "second order Godunov's schemes". I claim I cannot be considered as the 
author of them. Maybe, their authors were inspired by my suggestion about using the exact 
solutions of elementary problems to construct a numerical scheme. It is my pleasure to 
thank them for the attention they payed to me. 

Recently, two years ago, I learned that V. V. Ostapenko (Lavrentiev Institute of Hy- 
drodynamics SB RAS, Novosibirsk) studied actively the problem of an effective accuracy 
of numerical schemes in hydrodynamics based on asymptotic expansions. Instead of per- 
forming cumbersome analytical calculations, I advised him to make the experimental in- 
vestigations of the scheme we constructed together with V. S. Ryabenkii in 1957-1958. I 
hoped he would repeat our experiments and they would be described at last in details. But 
V. V. Ostapenko concentrated on the investigation of a more modern Lax-Harten scheme 
[19,21]. This scheme also demonstrated a mismatch between the formal order of approxi- 
mation and the effective order of accuracy (order of convergence). At the same time, from 
the experimental results of V. V. Ostapenko it apparently followed that the accuracy of the 
Lax-Harten scheme in the sense of weak convergence on discontinuous solutions is not less 
than the first order, i.e. the accuracy is higher than for my original scheme. 

I think, it is very interesting to carry out massive experimental investigations of pre- 
cision for all principal numerical methods. Especially interesting would be such a study 
for 2-D and 3-D problems which require the elaboration of the corresponding experimental 
techniques. 

One should remark the active work of K. A. Semendjaev to transform my scheme into 
a final computer code. His experience was very useful for our success. 
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3 2-D problems and moving grids 

Despite the fact that the work on the 1-D approach has been continued I started thinking on 
a 2-D variant of my algorithm. The attempts to solve gas dynamics problems were done in 
Steklov's Institute and later in our Institute by K. I. Babenko and I. M. Gelfand. Appear- 
ing difficulties were discussed on seminars where M. V. Keldysh took part. In particular 
during these works and accompanying discussions the numerical schemes - explicit in one 
spatial variable and impUcit in the other one ("sausages") - appeared. Also the first variant 
of the matrix double-sweep method that was suggested by M. V. Keldysh and investigated 
by K. I. Babenko and N. N. Chentsov appeared. There were a lot of discussions after the 
Locutsievskii's report based on the Rayleigh's book "Theory of Sound" chapter devoted to 
the instabiUty of a contact discontinuity. Indeed, the contact discontinuity evolves in time 
according to the law described by differential equations for which the Cauchy problem is 
not well posed in Hadamard's sense ( and not only its solutions are instable). The non well- 
posedness is characterized by the following fact - the short waves increase their amplitude 
very fast, the smaller the wave length is, the faster it is. The instability is a similar phe- 
nomenon, but in this case the characteristic time, during which the amphtude of the wave is 
increased, is bounded for all wave lengths. It was thought that processes described by the 
equations which are non well-posed in Hadamard's sense could not be computed numeri- 
cally because numerical schemes should be unstable. We discussed some regularizations 
of the boundary conditions on the contact discontinuity. However I cannot remember if any 
realistic variant of such a regularization was developed. 

Impressed by these discussions I got into a panic about the problems where contact 
discontinuities appeared. At that time there was a large interest in the problems that could 
be considered as near 1-D, i.e. the problems with slightly curved shocks and sound fronts 
and almost plane contact discontinuities. It occurred to me that for these problems a simple 
generaUzation of my scheme could be apphed, and I started developing such generaliza- 
tion. My solution was supported by M. V. Keldysh and I. M. Gelfand. K. A. Bagrinovskii 
and G. B. Alalykin took part in the development of the first variant of algorithm. Later, 
A. V. Zabrodin joined us. The codes were again developed by V. V. Lucikovich and G. 
N. Novozhilova. It is interesting to note that, when the 2-D computer code has been writ- 
ten, the problems became more complex due to increasing demands of engineering inter- 
ests. These complications led to modifications of the code and of the elements of the com- 
putational method. Busy with the computer code we forgot about the main danger - non 
well-posedness of the contact discontinuity. As far as I remember this non well-posedness 
did never appear. 

I cannot understand why we did not meet this non well-posedness. Apparently, finite- 
difference equations used for the calculation of the boundary trajectories provide some 
forced regularization. If it is really so, one should be anxious about the following. The 
mechanism of this regularization is unknown, an hence we are not able to guarantee that it 
does not produce some undesirable effects. 

From my point of view, the detailed analysis of the algorithms used for the calculation 
of contact discontinuities is still an actual problem. In this case, we can expect that some 
new effects will be discovered. 

One should note that calculations with moving grids, in which the tracked shock wave 
induces the displacement of points in 2-D and 3-D grids, has not been studied carefully, 
neither theoretically nor experimentally. The feedback of such point displacements has not 
been investigated. And it is not clear what kind of effects this feedback produces. 

When we implemented the 2-D approach based on the solutions of the Riemann prob- 
lem with arbitrary initial conditions, the first question concerned the construction of such 
solutions. In the 2-D case the rectangular grid cells can neighbor not only on each other 
but also on the nodes where four cells meet. If one constructs a 2-D scheme analogously to 
the 1-D, one should have analytic solutions of hydrodynamical equations with four discon- 
tinuities of initial data at one point. We did not have such solutions, even now they do not 
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exist, at least for general initial data. We had the audacity to suggest to use only classical 
solutions of the Riemann problem combined from the plane waves and describing initial 
Riemann problem placed on the edges of the neighboring cells. We ignored interaction of 
four cells having a common node. Implementing this approach we abandoned a clear phys- 
ical interpretation on which the construction of the 1-D scheme was based. Obviosly, there 
were a lot of discussions about the suggestion mentioned above. As a result we decided to 
perform the investigation with the help of a calculation of a shock or even an acoustic wave 
moving in the direction of a grid cell diagonal. At the same time, for the acoustic waves 
propagating in a medium at rest, K. V. Brushlinskii, based on Gelfand's suggestion, con- 
structed the solution of the problem using an interaction of all cells adjoint to one node with 
the help of the Sobolev's method of functionally invariant solutions. This solution was used 
in a numerical scheme completely analogous to the 1-D one. After that, the calculations of 
waves moving in a grid diagonal direction were performed with the help of Brushlinskii's 
scheme and compared with our scheme. To our surprise and satisfaction, we did not dis- 
cover any essential differences. After that, only the rough model was employed. We made 
a lot of efforts to derive a stability criterion needed to choose an admissible time step for a 
prescribed spatial step. First, using the Fourier method, G.B. Alalykin, K.A. Bagrinovskii 
and I, arrived at a cubic characteristic equation and derived from it only the necessary con- 
dition of stability. Then, together with K.A. Bagrinovskii, we succeeded in obtaining a 
sufficient condition of stabiUty considering results of 2-D calculation as some averaging of 

1- D calculations, or as one says by using a splitting approach. This work [24] was pub- 
lished in 1957 where a few not very simple tests were performed with the help of our 2-D 
code. 

The idea to use the splitting approach appeared under the influence of investigations of 

2- D implicit schemes for the heat equation. These schemes were analogous to ones that 
were proposed later by J. Douglas and their co-authors. But at that time we refused to use 
splitting for the 2-D heat equation. 

Similarly to the construction of the numerical schemes for the gas dynamics and acous- 
tics using non-smooth - discontinuous - solutions, I decided to test proposed spUtting 
schemes for the equation ut = Uxx + Uyy, taken as a test for the problem of the evolution 
of a heat impulse released in one grid cell. In other words, I decided to simulate the solution 
which after some time should not be essentially different from 



At least, the level curves of the computed solution should be convex curves - close to 
circles. 

I was very surprised when, during calculations with the Courant number 



I realized that these level curves turned out to be cross-like at least at the first time steps. 
This phenomenon forced me to avoid the spUtting for the heat equation. However, these 
investigations led us to use the splitting procedure not only for the stability analysis of gas 
dynamics numerical schemes but also for the organization of the code structure. 

The first 2-D code was constructed in such a way that on the successive steps the fluxes, 
computed with the help of the 1-D approach, were used at first in one direction and then 
in the other grid direction. There were no troubles, but later we refused to use splitting 
in these problems with exphcit schemes as V. V. Lucikovich suggested. He said that this 
refusal allowed to simphfy the code and made it more fast-acting. Before we were sure that 
using the splitting approach led to simplification. 

As I mentioned above we had to use the Lagrangian coordinates and turned to the Eule- 
rian coordinates. But at the same time in order to connect the grid with moving boundaries. 
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we had to make the moving grids as well. In particular, it allowed to include the shock 
fronts into a special type of boundary. The usefulness of this fact was demonstrated at the 
calculation of the flow around a sphere [6]. 

One should remark that in the first variants of our 2-D approaches we constructed 
meshes and wrote finite-difference formulas on the basis of conservation laws with a lot 
of caution and carefulness. For example, for cell boundaries, the arcs of logarithmic spi- 
rals were used, and the integrals over the cells limited by these spirals were calculated by 
explicit formulas. Only some years later we began using more simple variants. The hard 
work to derive the analytic formulas for integrals was performed by K. A. Bagrinovskii and 
A. V. Zabrodin. 

Even in 1 -D calculations that were performed in the Lagrangian coordinates we tried to 
increase the accuracy by looking at the advancement of the strongest shock waves. In order 
to avoid the spreading of these waves, their coordinates were marked and the grid cells 
where they were located were divided into two parts. To compute the displacement of a 
particular wave, the Hugoniot relations were used. After that, we started using the moving 
grids in the Eulerian coordinates, and the necessity to distinguish the boundary type - 
shock waves and the boundaries of material layers (contact discontinuity)- disappeared. 
It simplified the logical structure of the algorithm and allowed us to introduce one more 
boundary type on Zabrodin's suggestion. These boundaries had an assigned position or 
moved following a specified law and did not influence the medium flowing across them. 
They were called Eulerian boundaries. 

The grids in regions, Umited by boundaries of different types, were attached to these 
boundaries and computed with the help of simple interpolation formulas. Of cause, at each 
time step, the displacement of boundaries implied the displacement of grid points. The use 
of codes with such a multi-region structure allowed to increase the accuracy. It was possi- 
ble not due to the fact that numerical formulas were improved, but because the grid adapted 
to the solution's structure and strong discontinuities were not smeared. Working with this 
code I wanted to change the grid calculation based on the interpolating formulas in a ge- 
ometric region with assigned boundaries into the solution of some differential equations. 
These equations described a mapping transforming this region into some standard one (for 
example, into rectangular). G. P. Prokopov and 1 have studied this problem for seven years 
(1961-1968) to achieve the first acceptable results [5,25,26]. 

I am still interested in the classes of mappings for the grid generation although 1 stopped 
studying numerical schemes many years ago [27]. 

From the beginning, G. P. Prokopov and I postulated that the mapping should be defined 
as a solution of elliptic equation systems. Our main efforts were directed to solve them 
efficiently. We decided to use a variational approach that was realized with the help of a 
finite element method. However, we did not perform the first test on calculation of grids 
in hydrodynamical problems with moving boundaries. Because in this case the elliptic 
system should be solved on each time step (we were afraid of expensive time-cost). First, 
we decided to use a variational approach in some stationary problems and applied our ideas 
to calculate the critical parameters of a nuclear reactor (see [8]). I mentioned above that 
during these secondary investigations I was interested in computational methods of linear 
algebra. I was attracted by this problem, and I devoted to it my whole attention. 

Only after that we succeeded in finding the solution of stationary problems, we decided 
to generate grids on each step by solving elliptic equations. The first sufficiently universal 
code started to work at the end of 1968 - at the beginning of 1969. The first problem 
computed by it was based on the article [28]. This article contained the talk presented at 
the International conference on explosion physics (Novosibirsk, 1969) and was devoted 
to the wave formation at the explosion welding and the analysis of experiments that were 
carried out at the Novosibirsk Institute of Hydrodynamics. 

In September 1969, 1 moved to Novosibirsk and later did not devote myself to numerical 
schemes. 
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I should remark that during the development of the moving grids for 2-D calculations, 
we tried to use them in 1-D hydrodynamical problems to obtain more exact results. The 
whole domain was divided into some sub-domains whose boundaries were contact and 
shock discontinuities, and characteristics. Inside each sub-domain a scheme of the second 
order was used. It was important for the calculation of rarefaction waves to use many nodes 
immediately after they emerged (for instance, 50). It provided very high accuracy, and it 
was possible because of application of implicit numerical schemes. 

We carried out the described method together with 1. L. Kireeva and L. A. Pliner, and 
the algorithm and code were written by G. B. Alalikin. Based on this code some principal 
calculations were made. The method was published in 1970 (see book [1], in russian). 
I think the book passed unnoticed but 1 believe the speciahsts in the area of numerical 
hydrodynamics can find something interesting even today. 

One should mention that at the same time 1 worked with the characteristics method to 
adapt it to computer calculations. The result of this activity was the method of character- 
istics by layers in which not the intersection points of characteristics are calculated, but 
the coordinates of these characteristics at fixed layers t = const. But after having finished 
both the method and the code, during the first calculations the method showed its inade- 
quateness. It did not calculate strong rarefaction waves correctly. We did not expect such 
a situation because the characteristics method supposed to be good in calculation of such 
waves. Brief communication about this failure was published in [3] and it led me to con- 
struct a scheme of the second order to calculate solutions with distinguished discontinuities 
to which the above cited book is devoted [1]. 

The scheme for the gas dynamics calculations based on the solution of the Riemann 
problem began wide- spreading only in 1969, where even three independent presentations 
were reported during the conference in Novosibirsk. I have already mentioned our pre- 
sentation. Besides, there was the presentation of M. Ya. Ivanov and A. N. Kraiko from 
Moscow Central Institute of Aircraft Engines and the presentation of T. D. Taylor and V. 
S. Mason (USA) about using our scheme to calculate gas-dynamical flow around "Apollo" 
(see [29,30]). Later, together with our colleagues from Institute of Aircraft Engines we 
prepared a detailed description of the numerical method and published the book translated 
in French [31]. 

4 Conservation laws and thermodynamics 

I will touch one more issue appeared during the work on the construction of numerical 
schemes and related to the concept of generalized solutions of conservative quasilinear 
equations. 

As it was mentioned above the conservation laws of mass, momentum and energy are 
valid for the numerical solutions computed with the help of my scheme. However the 
classical solutions (without discontinuities) of gas dynamics equations obey one more ad- 
ditional conservation law - the entropy conservation law. The entropy increases in the case 
where gas passes through the shock front (discontinuity in the solution). This statement is 
the Zemplen's theorem, which represents (from the point of view of the theory of quasi- 
linear equations) the postulate included into the definition of generalized solutions. With 
the help of this postulate the discontinuous solutions satisfying the conservation laws of 
mass, momentum and energy, so that the entropy of a given material volume is decreasing, 
are excluded from the number of the generalized solutions. In my scheme, the law of non- 
decreasing of the entropy is automatically fulfilled because the solutions of the Riemann 
problem included in this scheme satisfy this law, and also the entropy is increasing in the 
process of averaging in computational cells. Such averages are performed at the end of 
each time step in my scheme. 

Obviously, 1 was interested in the description of such quasi-linear equations for n un- 
knowns whose classical solutions satisfy automatically to a n -|- 1-th conservation law. 
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Trying to understand which is the answer on this question, I noticed that all nonlinear rela- 
tions which are present in the hydrodynamical equations are defined by only one nonlinear 
function - the thermodynamic potential E{V, S). In the Gelfand's lectures on the theory 
of quasi-linear equations that he gave at the Moscow University [32], the role of in the 
equalities Es > and Evv Ess — Eyg > was mentioned. These inequalities in the 
thermodynamics are used for the gases and liquids if the parameters V, S are far from the 
phase-transition values. As a problem, I. M. Gelfand suggested to prove that these inequal- 
ities provide the well posedness of spreading of sound-waves by acoustic equations for 
a heat conducting gas. Impressed by this problem 1 supposed that similar considerations 
about well posedness or stability (for systems with finite numbers of degrees of freedom) 
can be put in the base of a well-known thermodynamic theorem on the existence of the in- 
tegrating multiplier 1/T for dE +pdV. This theorem is the basis of the entropy definition 
with the help of the equality dS = (l/T) (dE + p dV). After several months of hard work 
I succeeded to prove this hypothesis [33]. 1 found out that the stability, that proves the 
impossibility to construct the perpetual motion machine of the second kind, is always the 
consequence of existence of some integral which is dissipated due to the heat transfer. If 
we study only the stability of small oscillations (i.e. describing the evolution of a system by 
linear equations), the dissipating integral appears to be a quadratic form with a symmetric 
and positive definite coefficient matrix. The symmetry of this matrix leads to the equalities 
from which the existence of a universal integrating multiplier is derived. In my article [33] 
(see also § 22 in the book [23]) the drafts of the perpetual motion machine of the second 
kind are presented in the case where the universal integrating multiplier does not exist. 

After finishing the work [33], in order to understand the reasons for which the three 
equations have the "extra" fourth conservation law, it was natural to represent the hydrody- 
namical equations in the Lagrangian coordinates (they are used for classical solutions) 



dV_ 
'dt 



du 

dx 



= 



du ^ dEv{V,S) _^ 
dt dx ' 

and derive the fourth conservation law as their linear combination 



^^ d{E{V,S)+v?/2) ^ d{uEv{V,S)) 



= Ev{V,S) 



dt 

dV du 
dt dx 



+ u 



du 
dt 



+ 



dx 

dEv{V,S) 



dx 



+ Es 



95 
'dt 



Obviously, such an equality representing a linear dependence between four conservation 
laws can be solved relatively to any of them. In particularly, the entropy conservation law 

dS 

-— = can be presented as a linear combination of the conservation law for the specific 

dt 

volume 



dV du_^ 

dt dx 



the conservation law of the momentum 

du 



dEv _ Q 
dt dx ' 



and the conservation law of the energy in the following way 



Ey 

Es 



dV_ du 
dt dx 



u 

Es 



du dEv 

1 

dt dx 



1 

Es 



d{E + u^/2) d{uEv) 



dt 



+ 



dx 
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The last equality is possible due to the relation between the differentials 
dS = —^dV-—du+—diE+-). 

tjs -C'S -C/S ^ 



It is convenient to use the following form 
d 

and to introduce the notation 



Es 2Es Es 2 ^ 



V d-^ +ud— - {E + -)d— , 



L =S-±{E+'-() + ^ + ^ = ±[SEs+VEy-E], 

Ev u 11 

91 = ' 92 = t;- , 93 = = • 

£/S hig h/s 1 

As a result, we have the equalities that relate derivatives Lq. to the quantities appearing in 

d 

the hydrodynamical equations under the derivative — 

Lq,=V, Lq,=U, Lq,=-{E + U^/2) , 

and the relations between the differentials dS, dV, du, d{E + in an elegant form, 

on which the well-known Legendre's transform in the theory of convex functions is based 

d(9i W + 92 Lq^ + 93 Lq^ -L) = Lq^ dqx + Lq^ dq2 + Lq^ dgs . 
It is convenient to introduce one more function of 91, 92, 93 



such that 





M 


9192 






93 


Mq, 


- -Si - 


u , 


93 




Mq, 


- -Si - 


Ev , 


93 




Mq, 


9192 


= uEv , 


9i 




M - 


91 Mq, - 


92 Mq, - 



93 Mq,=0. 

In this case the hydrodynamical equations have the following form 

dLq dM„. 

They consist of three conservation laws, and the forth (additional) one 

d d 

— {qi Lq, + q2 Lq, + 53 Lq, -L)-— (gi Mq, + 52 Mg, + 53 M,3 - M) 

d 

= olil^ ^gi + 92 Lq, + q3Lq,-L) = 0, 

dS 

coincides with the conservation law of entropy = 0. Afterward it was not difficult to 

obtain for other classical equations of mathematical physics, even multidimensional, the 
standard representation 

«1 «L = . 

dt dxj 
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I succeeded in writing in a similar form the gas dynamics equations in the Eulerian co- 
ordinates. I was happy with the fact that if the generating "thermodynamic potential" 
L = L{qi,q2, . . .) is convex, theseequationscanberewrittenin theformof theFriedrichs's 
symmetric hyperbolic system 

dqk dqk _ 

with a positive definite matrix of coefficients at derivatives dqt on t. K. O. Friedrichs 
showed that for such systems the Cauchy problem is well-posed if some conditions on the 
smoothness of initial data are true. 

So, my hopes to discover the connection between well-posedness and the laws of phe- 
nomenological thermodynamics have been realized. 

I described all my achievements in the article that was sent to the journal "Uspekhi 
Matematicheskikh Nauk", but it was not published because a reviewer decided that it did 
not have a mathematical content. 

After that 1 was dealt a new severe blow. 1 decided to use the obtained equations by 
introducing into them small additional dissipation terms. 1 wanted to confirm that the dis- 
continuous solutions, after introducing any small viscosities, are just "smeared" in a thin 
stripe, and that the limit solutions (generalized solutions of equations with zero viscos- 
ity) do not depend on the form of these viscosities. Considering solutions in the form of 
travelling waves qt = qi{^),^ = {x — u)t)/e for equations 



dig, _^ dMg^ _ d 
dt dx dx 



1 could propose an obvious geometric interpretation for ordinary differential equations 
for qi{£,). This interpretation led to a construction of systems of hyperbolic conservation 
laws for which admissible discontinuous solutions satisfied the entropy increasing law, and 
which depend on what kind of small dissipative terms we introduced to smooth them. In 
physical problems the real dissipative processes can differ from the numerical dissipation 
that depends on the numerical scheme. Soon, B. F. Diachenko constructed the example in 
which different numerical approximations led to different numerical solutions [34]. 

1 was very upset that 1 could not justify all hypotheses lying in the background of 
the numerical scheme. My work with the description of this example was presented in 
"Doklady AN SSSR" [35], [36] by 1. G. Petrovskii. Before the article was published, 
1 presented it on the seminar where R. Courant and R D. Lax participated, and visited 
Moscow for the first time. 

Later, in Novosibirsk, together with my colleagues 1 continued the investigations on 
conservation laws and their connections with the thermodynamics. The review of a part of 
these investigations was included into the paper [22] that I presented in 1986 in St.-Etierme 
at the Conference on hyperbolic equations. 1 should also mention important studies [37], 
[38] on this issue performed by R D. Lax and K. O. Friedrichs. Last years these studies 
were continued and presented by me and E.l.Romenskii in Lisbon, Lake Tahoe (USA, 
Nevada), Paris [39-41]. Our work [42] (in collaboration with T. Yu. Mikhailova) was 
also devoted to this topic. The works [41], [42] showed unexpected connection between 
formally overdetermined systems of conservation laws in mathematical physics and the 
theory of representations of the rotation group. 

1 should notice a recent work [43] in which the estimates of the entropy growth for the 
systems of conservation laws in the form described above are investigated. The author of 
this work could connect these estimates with the variation of solutions. It seems to me that 
this original way of definition of such a variation is very promising. For this definition it is 
not important if the problem is one-dimensional or multi-dimensional. 
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5 Conclusion 

In this note I tried to describe the intensive work in the field of numerical hydrodynamics 
in which 1 was involved by my supervisers I. G. Petrovskii and I. M. Gelfand, and the result 
was the development of Godunov's scheme. 

I should emphasize the influence of the whole scientific community within I worked. 
Of course, different approaches to the numerical hydrodynamics were elaborated by many 
other research groups and some of these groups belonged to the Institute where I worked. 
The scientists who participated in these groups are K. I. Babenko, V. V. Rusanov, I. M.Gel- 
fand, V. F. Diachenko, A. A. Samarskii, N. N. Janenko, B. L. Rozhdestvenskii, and others. 
My goal was not to describe the whole history of this subject but only present the part in 
which I participated. 

Various discussions with outstanding physics theoreticians of that time such as Ya. B. 
Zeldovich, A. D. Sakharov, D. A. Frank- Kamenetskii, Yu. B. Hariton, as well as members 
of the experimental group of L. V. Altshuler were very essential for me. I have already 
remarked that the comments of physicists were, as a rule, very critical. They stimulated the 
detailed analysis of all numerical effects and their reasons, and led to modifications of the 
method. The specific problem for a steady flow around a body was suggested to me by G. 
I. Petrov. 

I wanted to demonstrate how the set of the principal scientific issues emerged during 
that initial period of the development of the computational fluid dynamics. Many of these 
issues are stiU unresolved even today. 
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Appendix 

Experimental study of a disparity between order of approximation 
and accuracy 

We present a method allowing an efficient estimation of the order of the weak con- 
vergence, which is based on generalized solutions of hyperbolic systems of conservation 
laws. The idea of the method was suggested by S. K. Godunov and V. S. Ryabenkii. The 
method is based on the experimental estimation of the convergence rate of the first inte- 
grals of numerical solutions which are calculated over domains with singularities of the 
approximated exact solution. The analysis of the accuracy of the Harten's TVD-scheme 
with the second order approximation on smooth solutions shows that only the first order of 
convergence occurs in the problem of dam breaking with creation of a bore and a smooth 
rarefaction wave. It results in the decreasing of the order of strong convergence in smooth 
parts of exact solutions behind the wave front. The result obtained is in contradiction with a 
wide-spread opinion according to which the monotonic numerical schemes of higher order 
of approximation on the smooth solutions conserve a higher order of local convergence on 
smooth solutions of quasi-linear hyperbolic conservation laws. 

Let us consider the Cauchy problem for hyperbolic system of quasi-linear conservation 
laws 

ut + f{u)x^O, u{0,x) = uo{x) , xeR, (1) 

where uo{x) is a m-dimensional piecewise continuous vector-function and f{u) a smooth 
flux function. We suppose that the problem (1) has a unique stable generalized solution 
u{t, x). The numerical solution is calculated with an explicit two-layer in time and sym- 
metric in space conservative scheme 

= ^" - ^" (/JVi/2 - /r-1/2) , v'^^uoijh), jeZ, (2) 

where /J^./^ - . . . , 1;!^+ J, f{u,...,u) = f{u), for all u 6 W'\ = 

v{tn,j h), to = 0, tn = To + Ti + ■ ■ ■ + r„_i, A„ — Tn/h, Tn is the time step at n-th time 
layer t„, h is the constant spatial step, / is the continuous function of numerical flux. The 
time step t„ is calculated from the Courant condition r„ < r™^'' = z h/ maxij |a*(w")|, 
where z — 0.5 is a safety factor, a'(u), i — 1, . . . ,in are the eigenvalues of Jacobi matrix 
fu in (1). Let us fix a number a, a G R and introduce the integrals 

U'^iT.x)^ u{T,y)dy, V,:{T,x)^ Vh{T,y)Ay . (3) 

J a J a 

The numerical solution Vh{T,x) converges weakly to the exact solution u{T,x) with r- 
order within segment [a,x] C K if 

V^iT, x) - [/"(T, x) = Ch'- + o{h'), (4) 

where C does not depend on h. 

If [a, x] does not contain singularities of the exact solution, the order of weak conver- 
gence coincides with the order of a local convergence on smooth solutions, i.e. equals to 
2 in our case. Otherwise, the order of weak convergence is less than 2, because the TVD 
schemes do not have second order of weak approximation on discontinuous functions. 

Let us estimate the order of weak convergence r in the case where an exact discontin- 
uous solution of the problem (1) is unknown in advance. To calculate the order of weak 
convergence (based on the Runge's rule) it is enough to have three numerical results with 



'In this Appendix the results obtained by V.V. Ostapenko (Lavrentev Institute of Hydrodynamics SO RAS, 
Novosibirsk) are used. 
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rather small spatial steps hi = h,h2 = h/2, h^ = h/A. For each step the property (4) is 
satisfied, hence 

sv^ = IK - I - \c\ (K K+i) + o{hn .. * - 1, 2. (5) 

Therefore, we obtain with the accuracy o{h^) 

5Vi _h\-hl _ 1 - (1/2)'- 



5V2 hi - hi {\/2Y ~ (1/4)' 
from which it follows that 



1 -^^i (a^ 
r = lo&^. (6) 

Below we present the results of a numerical calculation of errors (5) and the order of 
accuracy (6) based on the TVD-scheme in the problem of dam breaking in a channel of 
non-constant depth. 

As a test scheme one chooses one of the TVD-schemes suggested by Harten 0. This 
scheme is applied to the numerical solution of the Saint- Venant equations (shallow water 
equations). They are equivalent to the equations of the polytropic gas with polytropic 
exponent 7 = 2 and have the following form 

H \ f q 



where H and q are the depth and the flux, respectively, g is the gravity (in calculations 
9 = 10). 

In the Figures 1-4 for T = 0.2 one can observe the results of the calculation of the 
Cauchy problem with discontinuous initial data 

^ ' ^ I 2 - th{x - 5) , X > 3.5 , ^ ^ 

q{0,x)^0, xeR. (8) 

In these calculations, hi = h = 0.1, /i2 = h/2 = 0.05, h^ = h/A = 0.025. 

The results of the calculations are presented in Figure 1 : circles correspond to the basic 
grid with step hi, the solid line to the grid with step h^; the dotted line corresponds to the 
initial position of the water level given by function (7). 

Figures 2-3 are the error, as defined by (5), and the order of accuracy given by relation 
(6); a = 12 in relation (3). From these plots, one can see that, when x > Q, i.e. when the 
segment of integration [a, x] lies entirely in the smooth part of the solution just before the 
shock wave front, the order of weak convergence r{x) equals 2. When x < 5 and when 
[a, x] contains the front of a discontinuous wave, then r{x) « 1, and the TVD-scheme (2) 
has approximately only the first order of accuracy for the Cauchy problem (7)-(8). 

In Figure 4, the function 



is plotted where 



R{x) , R{x) < 3 
3 , R{x) > 3 



R{x) = log ^"^1^1 , 6v,{x)^\vh^{T,x)-Vh^ {T,x)\, i = l,2, 

0V2{X) 

that presents the order of the local strong convergence of the numerical solution. As one can 
see from Figure 4, the calculation domain is divided into three isolated parts by singularities 
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of the exact solution (inside every part the order of local convergence is changed in a 
random way). Inside each part there is a different order of the local convergence. The order 
of the local convergence of the TVD-scheme on different smooth parts is quite different, 
and, generally speaking, less than the order of approximation on the smooth solutions (see 
Fig.4). 




Figure 1: Comparison of an exact and a numerical solution 
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Figure 3: The weak convergence order r 
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Figure 4: The strong convergence order ro (x) 
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